function [ v, zeta , Lh , L , vL , C ] = ...
    SolveC( Grids , Params , GE )

% This function solves for C using bisection

ScaleFactor = Grids.ScaleFactorC ;
bargaining   = Params.barg ;
x0 = Params.muLearn / ( Params.muLearn + Params.DeltaLearn ) ;                

% Define function that updates C
run MiscFunctions.m
  
% Initialize bounds

C0 = GE.C ;
C1 = GE.C ;

% Lower bound
GE.C = C0 ;
[ v , zeta , Lh , L , vL ] = SolveVL( Grids , Params , GE ) ;
Cratio0 = EntryRatio(x0,v,zeta,Grids,bargaining,GE.C) ;
enteredloop0 = 0 ;
% Cratio is a decreasing function of C

if abs(Cratio0 - 1 ) > Grids.tolC 

    while Cratio0 < 1
        C0 = C0 / ScaleFactor ;
        GE.C = C0 ;
        [ v , zeta , Lh , ~ , ~ ] = SolveVL( Grids , Params , GE ) ;
        Cratio0 = EntryRatio(x0,v,zeta,Grids,bargaining,GE.C) ; 
        enteredloop0 = 1 ;
    end

    % Upper bound
    if enteredloop0 == 1
        C1 = ScaleFactor * C0 ;
        GE.C = C1 ;
    else
        GE.C = C1 ;
    end
    [ v , zeta , Lh , ~ , ~ ] = SolveVL( Grids , Params , GE ) ;
    Cratio1 = EntryRatio(x0,v,zeta,Grids,bargaining,GE.C) ; 

    enteredloop1 = 0 ;      
    while Cratio1 > 1
        C1 = ScaleFactor * C1 ;
        GE.C = C1 ;
        [ v , zeta , Lh , ~ , ~ ] = SolveVL( Grids , Params , GE ) ;
        Cratio1 =  EntryRatio(x0,v,zeta,Grids,bargaining,GE.C) ;
       enteredloop1 = 1 ;
    end

    % bisection loop
    if enteredloop1 == 1
        C0 = C1 / ScaleFactor  ;
    end

    error = 1 ;
    it = 1 ; 
    while error > Grids.tolC && it < Grids.maxit

        Cmid = ( C0 + C1 ) / 2 ;
        GE.C = Cmid ;
        [ v , zeta , Lh , ~ , ~ ] = SolveVL( Grids , Params , GE ) ;
        Cratio = EntryRatio(x0,v,zeta,Grids,bargaining,GE.C) ;

        error = abs( Cratio - 1 ) ;
        disp(['   Bisecting on C: Cratio = ' num2str(Cratio) ])

        if Cratio < 1
            C1 = Cmid ;
            Cratio1 = Cratio ;
        elseif Cratio > 1
            C0 = Cmid ;
            Cratio0 = Cratio ;
        end

        it = it + 1 ;

    end

    frac = ( Cratio - Cratio0 ) / ( Cratio1 - Cratio0 ) ;
    C = ( 1 - frac ) * C0 + frac * C1 ;
    GE.C = C ;
    [ v , zeta , Lh , vL , ~ ] = SolveVL( Grids , Params , GE ) ;

else
    
    C = GE.C ;

end

% Recover actual population density
L = Lh / sum( Lh .* Grids.mx ) ;


end

